!=======Programa que calcula movimento do pendulo por Euler-Cromer======

PROGRAM exerA2

!=======================================================================

IMPLICIT NONE

!=======================Declaração de variáveis=========================

      INTEGER :: I, steps
      REAL(8) :: tethai, tethai1, omegai, omegai1, Dt, currt
      REAL(8) :: m, l, totalt, emec
      REAL(8), PARAMETER :: g = 9.80665d0, PI = 3.14159265358979d0

!=============Inicialização, leitura e abertura de arquivos=============

      WRITE(*,*) "Entre com a massa na ponta do pêndulo:"
      READ(*,*) m
      WRITE(*,*) "Entre com o comprimento da haste:"
      READ(*,*) l
      WRITE(*,*) "Entre com o intervalo entre iterações delta_t:"
      READ(*,*) Dt
      WRITE(*,*) "Entre com o ângulo inicial do pêndulo (em graus):"
      READ(*,*) tethai
      WRITE(*,*) "Entre com o tempo total de simulação:"
      READ(*,*) totalt

      !Abre arquivo para tetha(t)
      OPEN(UNIT=10, FILE='exerA2_out.dat', STATUS='unknown')
      !Abre arquivo para E(t)
      OPEN(UNIT=20, FILE='enerA2.dat', STATUS='unknown')

      !Converte o angulo de graus para radianos
      tethai = tethai*PI/180.0d0

      tethai1 = 0.0d0
      omegai = 0.0d0
      omegai1 = 0.0d0
      currt = 0.0d0

      !Calcula a energia inicial do pêndulo (só potencial)
      emec = 0.5d0*m*g*l*(tethai**2.0d0)

!=============================Faz as iterações==========================

      steps = totalt/Dt
      DO I=1,steps
            !Escreve no arquivo
            WRITE(10,*) currt, tethai
            WRITE(20,*) currt, emec
            !Faz a iteração de omega
            omegai1 = omegai - (g/l)*tethai*Dt
            !Faz a iteração de tetha
            tethai1 = tethai + omegai1*Dt
            !Testa para ver se tetha está nos limites
            IF (tethai1 .GE. PI) THEN
                  tethai1 = tethai1 - 2.0d0*PI
            ELSE IF (tethai1 .LE. -1.0d0*PI) THEN
                  tethai1 = tethai1 + 2.0d0*PI
            END IF
            !Calcula a energia
            emec = 0.5d0*m*(omegai1**2.0d0)*(l**2.0d0) +&
            &0.5d0*m*g*l*(tethai1**2.0d0)
            !Calcula em qual tempo ocorreu essa iteração
            currt = currt + Dt
            !Prepara variavaies para a proxima iteração
            omegai = omegai1
            tethai = tethai1
      END DO

!=======================================================================

END PROGRAM exerA2

!=======================================================================
